

 %
 % Related parameters:
omega0  = 1-beta;
omega1  = beta*alpha/(1 - beta*(1-alpha));  
mu      = -.5*(sigma^2)*alpha/(1-alpha);
mulog   = mu - ((sigma^2)/2);



%% NATURAL WASTAGE
muhat= mu + (1-alpha)*s*lambda;
rho0 = r + s*lambda;
rho1 = rho0 - muhat;
bb   = muhat - .5*(sigma^2);
cc   = 2*(sigma^2)*rho0;
gamma1 = (-bb - sqrt(bb^2 + cc)) / (sigma^2);
gamma2 = (-bb + sqrt(bb^2 + cc)) / (sigma^2);
G0 = 2.5;	
foptions = optimoptions('fsolve','Display','off');
[G,fval,exitflag] = fsolve('solve_nwfun', G0, foptions,  c,rho0,omega0,gamma1,gamma2);
if exitflag <1
    disp('WARNING: Failed to solve natural wastage boundaries');
    pause;
else
    thetaG = (G.^gamma2  - G)  ./ (G.^gamma2  - G.^gamma1 );
    phiG   = 1 - (thetaG /gamma1) - ((1-thetaG) /gamma2);
    mL = (omega0/rho0) / (phiG*(1-omega1)/rho1);
    mM = mL*G;
    J1 =    -thetaG *((1-omega1)/rho1)*((mL^(1-gamma1))/gamma1);
    J2 = -(1-thetaG)*((1-omega1)/rho1)*((mL^(1-gamma2))/gamma2);
end



%% UPPER SUPPORT IN NET GROWTH REGION, delta(mU) =0
delta_mM   = s*lambda;
const      = ((1-omega1)*mM/alpha) - omega0 - (r+delta_mM)*c;
[mU,fUval, exitflag] = fzero(@(x) (s*lambda + (1/c)*( (((1-omega1)*(x - mM))/alpha) - const*( (x/mM).^(1/(1-alpha)) - 1 ) )), mM*1.1 ); 
if exitflag <1
    disp('ERROR: Failed to solve for mU');
    pause;
end
if mU <mM
    disp('ERROR: mU < mM');
    pause;
end


%% DENSITY of employees over m
solve_ssdensity;


%% POLICY FUNCTIONS
solve_policy;


%% SS DISTRIBUTION of Marg. Product over **Employers**
solve_ssdensity_estab;



    
    
    